How to plot Pressure and Wind¶
- Run the cell below before running any other cell.Cell is commented. Uncomment it first before running it.
- Then restart the kernel and run other cell
- Keep in mind to perform these procedures each time you open this notebook.
In [1]:
# !apt-get install -qq libgdal-dev libproj-dev
# #!pip install --no-binary shapely shapely
# !pip install --no-binary shapely shapely --force
# !pip install cartopy
In [2]:
import xarray as xr
from datetime import datetime, timedelta
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings("ignore")
In [3]:
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
About Data Download¶
Now we will use data from "ERA5 hourly data on single levels from 1940 to present" that is availble for free in climate data store website: cds.climate.copernicus.eu
The data I downloaded has the following specifications:
Product type:
Reanalysis
Variable:
100m_u_component_of_wind,
100m_v_component_of_wind,
10m_u_component_of_wind,
10m_v_component_of_wind,
surface_pressure
Year:2023
Month:June, July
Day:
01, 02, 03, 04, 05, 06, 07, 08, 09, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31
Time:
00:00, 01:00, 02:00, 03:00, 04:00, 05:00, 06:00, 07:00, 08:00, 09:00, 10:00, 11:00, 12:00, 13:00, 14:00, 15:00, 16:00, 17:00, 18:00, 19:00, 20:00, 21:00, 22:00, 23:00
Sub-region extraction:
North 29°, West 85°, South 18°, East 97°
Format:
NetCDF (experimental)We can download the data directly using API. we need to follow the setps below to download the data:
Install the CDS API: Run the snippet below¶
In [4]:
!pip install cdsapi
Run the the snippet below to download the data:¶
Don't forget to change the output_folder location as you desire.
In [5]:
import cdsapi
import os
c = cdsapi.Client()
output_folder= '/mnt/102ECBA62ECB8368/Academic journey/ML_HEATWAVE/Data/Test/'
#Ensure theoutput folder exists
os.makedirs(output_folder,exist_ok=True)
c.retrieve(
'reanalysis-era5-single-levels',
{
'product_type': 'reanalysis',
'format': 'netcdf',
'variable': [
'100m_u_component_of_wind', '100m_v_component_of_wind', '10m_u_component_of_wind',
'10m_v_component_of_wind', 'surface_pressure',
],
'year': '2023',
'month': [
'06', '07',
],
'day': [
'01', '02', '03',
'04', '05', '06',
'07', '08', '09',
'10', '11', '12',
'13', '14', '15',
'16', '17', '18',
'19', '20', '21',
'22', '23', '24',
'25', '26', '27',
'28', '29', '30',
'31',
],
'time': [
'00:00', '01:00', '02:00',
'03:00', '04:00', '05:00',
'06:00', '07:00', '08:00',
'09:00', '10:00', '11:00',
'12:00', '13:00', '14:00',
'15:00', '16:00', '17:00',
'18:00', '19:00', '20:00',
'21:00', '22:00', '23:00',
],
'area': [
29, 85, 18,
97,
],
},
os.path.join(output_folder,'sp_u_v.nc')
)
In [6]:
data_path='/mnt/102ECBA62ECB8368/Academic journey/ML_HEATWAVE/Data/Test/sp_u_v.nc'
ds=xr.open_dataset(data_path)
In [7]:
ds
Out[7]:
expver = 1 means ERA5 and expver = 5 means ERA5T.
Google to understand more.
Check the dimensions¶
In [8]:
ds.sp.dims
Out[8]:
Here we will use just era5(expver=1) data and we will use data for time = 2023-06-15T12:00:00
In [9]:
expver_slice=1
desired_time = '2023-06-15T12:00:00'
surface_pressure = ds.sel(expver=expver_slice,time=desired_time)
In [10]:
surface_pressure
Out[10]:
Now we have data for 12PM,June 15,2023. We can plot this data in two ways
- Directly from the data
- Creating a dataframe and plotting
Plotting Directly from data¶
In [11]:
from matplotlib.cm import get_cmap
warnings.filterwarnings("ignore")
# Get latitude and longitude from xarray.Dataset
latitude = surface_pressure['latitude'].values
longitude = surface_pressure['longitude'].values
# Get sp,u10 and v10 from xarray.Dataset
pressure_data = surface_pressure.sp
u10_data = surface_pressure.u10
v10_data = surface_pressure.v10
# Calculate wind speed and direction
wind_speed = np.sqrt(u10_data**2 + v10_data**2)
wind_direction = np.arctan2(v10_data, u10_data) * (180 / np.pi)
# Create a Cartopy PlateCarree projection
projection = ccrs.PlateCarree()
# Create a figure and axis
fig, ax = plt.subplots(subplot_kw={'projection': projection}, figsize=(12, 6))
# Plot surface pressure using contourf
contourf = ax.contourf(longitude, latitude, pressure_data, levels=10, cmap=get_cmap("jet"),transform=projection)
# Plot isobar lines using contour
contour = ax.contour(longitude, latitude, pressure_data, levels=10, colors='black', transform=projection)
# Plot wind direction using quiver
quiver = ax.quiver(longitude, latitude, u10_data, v10_data, wind_direction, transform=projection, pivot='middle', color='black')
# Add colorbar
cbar = plt.colorbar(contourf, ax=ax, orientation='vertical', pad=0.08, label='Surface Pressure (Pa)')
# Set title and labels
ax.set_title('Surface pressure and Wind Direction')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
# Add coastlines and gridlines
ax.coastlines()
ax.add_feature(cfeature.BORDERS)
ax.gridlines(draw_labels=True)
# Show the plot
plt.show()
Creating a dataframe and plotting¶
In [12]:
import pandas as pd
import numpy as np
sp=np.array(surface_pressure.sp)
u10=np.array(surface_pressure.u10)
v10=np.array(surface_pressure.v10)
latitude=surface_pressure['latitude'].values
longitude=surface_pressure['longitude'].values
reshaped_sp=sp.reshape(-1)
reshaped_u10=u10.reshape(-1)
reshaped_v10=v10.reshape(-1)
df = pd.DataFrame({
'latitude': np.repeat(latitude, len(longitude)),
'longitude': np.tile(longitude, len(latitude)),
'sp': reshaped_sp,
'u10': reshaped_u10,
'v10': reshaped_v10
})
In [13]:
df
Out[13]:
In [14]:
# Get the latitude and longitude values
latitude=df['latitude'].unique()
longitude=df['longitude'].unique()
# Reshape the sp,u10,v10 data into a 2D array matching the grid
sp_array=df['sp'].values.reshape(len(latitude),len(longitude))
u10_array=df['u10'].values.reshape(len(latitude),len(longitude))
v10_array=df['v10'].values.reshape(len(latitude),len(longitude))
# Create a Cartopy PlateCarree projection
projection = ccrs.PlateCarree()
# Create a figure and axis
fig, ax = plt.subplots(subplot_kw={'projection': projection}, figsize=(12, 6))
# Plot wind speed using contourf
contourf = ax.contourf(longitude, latitude, sp_array, levels=10, cmap=get_cmap("jet"),transform=projection)
# Plot isobar lines using contour
contour = ax.contour(longitude, latitude, sp_array, levels=10, colors='black', transform=projection)
# Plot wind direction using quiver
quiver = ax.quiver(longitude, latitude, u10_array, v10_array, wind_direction, transform=projection, pivot='middle', color='black')
# Add colorbar
cbar = plt.colorbar(contourf, ax=ax, orientation='vertical', pad=0.08, label='Surface Pressure (Pa)')
# Set title and labels
ax.set_title('Surface pressure and Wind Direction')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
# Add coastlines and gridlines
ax.coastlines()
ax.add_feature(cfeature.BORDERS)
ax.gridlines(draw_labels=True)
# Show the plot
plt.show()
In [ ]: